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Abstract 

This paper deals with the optimization of trial states for the computation 
of dominant eigenvalues of operators and very large matrices. In addition 
to preliminary results for the energy spectrum of van der Waals clusters, 
we review results of the application of this method to the computation of 
relaxation times of independent relaxation modes at the Ising critical point 
in two dimensions. 

I. INTRODUCTION 

The computation of eigenvalues and eigenstates of operators and large matrices is a 
ubiquitous problem. In this paper we review recent applications of the Quantum Monte 
Carlo methods that we have developed for this purpose. The reader is referred to other 
papers for introductory or more technical discussions of earlier work. [|T]-^| 

II. MATHEMATICAL PRELIMINARIES 

For an operator G, the power method can be used to compute the dominant eigenstate 
and eigenvalue, \ip ) and A . 

This well-know procedure can be summarized as follows: 

1. Choose a generic initial state \u^) of the appropriate symmetry. 

2. Iterate: 

\u^) = —G\u®), (1) 

where c t puts \u®) in standard form. 
For projection time t — > oo the following is almost always true: 
1. Eigenstate: 

M t+1) > - |V>o> (2) 
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2. Eigenvalue: 



A, 
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To see this, expand the initial state in normalized eigenstates 



u<°>> = X>?>hfc) 



(4) 
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with spectral weights wf\ Then \vfi) has spectral weights 
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This method can be implemented by means of a Monte Carlo method and, unlike varia- 
tional Monte Carlo, it has the advantage of producing unbiased results for large projection 
times t. The disadvantage is, however, that at the same time the statistical noise increases 
exponentially, unless G is a Markov (stochastic) matrix, or can be explicitly transformed to 
one. The statistical errors grow with the extent to which G fails to conserve probility, and 
to alleviate this problem, approximate dominant eigenstates can be used. 

In the case of Markov matrices, computation of the dominant eigenvalue is of no interest, 
since it is equals unity, but sampling the corresponding eigenstate has numerous applications. 



Given a set of basis states, one can construct trial states as linear combinations to 
obtain approximate excited or, more generally, sub-dominant states and the corresponding 
eigenvalues. These are computed by solving a linear variational problem. In a Monte 
Carlo context, the Metropolis method can be used to evaluate the required matrix elements. 
Subsequently, a variation of the power method can again be used to remove systematically 
the variational bias. [f)|-[7j Again, the price to be paid for reduction of the variational bias 
is increased statistical noise, a problem that can be mitigated by the use of optimized trial 
states. 

The linear variational problem to be solved for the computation of excited states is the 
following one. Given n basis functions find the n x n matrix of coefficients D^' such 
that 



are the "best" variational approximations for the n lowest eigenstates of some Hamilto- 
nian TC. In this problem we shall, at least initially, use the language of quantum mechanical 
systems, where one has to distinguish the Hamiltonian from the imaginary-time evolution 
operator G = exp(— tTC). In the statistical mechanical application discussed below, we 
shall encounter only the equivalent of the latter, which is the Markov matrix governing the 
stochastic dynamics. In the expressions to be derived below, the substitution 7iG p — > G p+1 



A. Subspace iteration 



n 



(6) 
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will produce the expressions required for the statistical mechanical application, at least if 
we assume that the non-symmetric Markov matrix that appear in that context has been 
symmetrized, which can always be accomplished if detailed balance is satisfied. 

Given these basis states, one seeks the "best" solution to the linear variational problem 
in Eq. (0) in the sense that for all i the Rayleigh quotient (ipi\T-C\tpi) / is stationary 

with respect to variation of the coefficients of the matrix D. The solution is that the matrix 

(i) 

of coefficients D\ has to satisfy the following generalized eigenvalue equation 

n n 

H ki D { P = Ej'ENuDjp, (7) 
i=i i=i 

where 

H k i = (u k \H\ui), and N ki = (u k \ui). (8) 

We note a number of important properties of this scheme. Firstly, the basis states |tt.j) in 
general are not orthonormal. Secondly, it is clear that any nonsingular linear combination of 
the basis vectors will produce precisely the same results, obtained from the correspondingly 
transformed version of Eq. (p]). The final comment is that the variational eigenvalues bound 
the exact eigenvalues from above, i.e., Ei > E iy where we assume E\ < Ei < . . .. One 
recovers exact eigenvalues Ei and the corresponding eigenstates, if the \ui) span the same 
space as the exact eigenstates, or in other words, have no admixtures of more than n states. 

The required matrix elements can be computed using the standard variational Monte 
Carlo method. The power method can subsequently be used to reduce the variational bias. 
Formally, one simply defines new basis states 

\u[ p) ) = G p \ Ui ) (9) 

and substitutes these new basis states for the original ones. In quantum mechanical appli- 
cations, where G = exp(— tH), the corresponding matrices 

h$ = (^\n\u?) (io) 

and 

= (4 P V/ P) ) (ii) 

can be computed by pure-diffusion Monte Carlo. || We note that, Monte Carlo yields these 
matrix elements up to an irrelevant overall normalization constant. 

As an explicit example illustrating the nature of the Monte Carlo time-averages that 
one has to evaluate in this approach, we write down the expression for N-j as used for the 
computation of eigenvalues of the Markov matrix relevant to the problem of critical slowing 
down, discussed in detail in the next section. One estimates this matrix as 

where the S t are configurations forming a time series that is designed to sample the distri- 
bution of a system in thermodynamic equilibrium, i.e., the Boltzmann distribution ip^. It 
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turns out that in this particular case, this distribution, the dominant eigenstate, has suffi- 
cient overlap with the magnitude of the sub-dominant states so that one can compute all 
matrix elements iV^ simultaneously without introducing a separate guiding function 0. 

The expression given in Eq. fll2|) yields the w/^B-auto-correlation function at lag p. 
The expression for H-?' is similar, and represents a cross-correlation function involving the 
configurational eigenvalues of the Markov matrix in the various basis states. Compared 
to the expressions one usually encounters in applications to quantum mechanical problems, 
Eq. ( |i~2D takes a particularly simple form in which products of fluctuating weights are absent, 
because one is dealing with a probability conserving evolution operator from the outset in 
this particular problem. 

III. UNIVERSAL AMPLITUDE RATIOS IN CRITICAL DYNAMICS 

Before continuing our general discussion, we temporarily change the topic to introduce 
stochastic dynamics of critical systems. What make such systems interesting, is that one 
can distinguish universality classes in which behavior does not depend on many of the 
microscopic details. For static critical phenomena, it is known that universality classes can 
be identified by dimensionality, symmetry of the order parameter, and the range of the 
interactions. For dynamical phenomena, there are additional features such as whether or 
not the dynamics is local or subject to conservation laws. 

On approach of a critical point, the correlation length £ diverges. The dynamical ex- 
ponent z governs the corresponding divergence of the correlation time r by means of the 
relation t a f . Since the critical exponent z is one of the universal quantities, it has 
been used to identify universality classes. Unfortunately, z does not vary by much from one 
universality class to another, and this poses a serious computational problem in terms of 
the accuracy required to obtain significant differences. One of the outcomes of the work 
reviewed here is that there are other quantities within computational reach, namely univer- 
sal amplitude ratios. M These ratios may serve as additional, and possibly more sensitive 
identifiers of universality classes. We shall consider various systems belonging to a single 
universality class, and we assume that the representatives of the class are parameterized by 

K. 

If a thermodynamic system is perturbed out of equilibrium, different thermodynamic 
quantities relax back at a different rates. More generally, there are infinitely many indepen- 
dent relaxation modes for a system in the thermodynamic limit. The Monte Carlo methods 
reviewed here have been used to compute relaxation times of Ising models on square L x L 
lattices at the critical point. || 

Let us denote by T K i(L) the relaxation time of mode i of a system of linear dimension 
L. As indeed scaling theory suggests, it turns out that the relaxation time has the following 
factorization property 

r Ki {L) « m K AiL z , (13) 

where m K is a non-universal metric factor, which differs for different representatives of the 
same universality class as indicated; is a universal amplitude which depends on the mode 
i; and z is the universal dynamical exponent introduced above. 
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Formulated as a computational problem, one has the following. Suppose S = 
(si, sip), with Si = ±1, is a spin configuration and p*(S) is the probability of finding S 
at time t. The probability distribution evolves in time according to 

Pm(S) = £P(S|S')MS'). (14) 

S' 

The detailed structure of the Markov matrix P is of no immediate importance for the 
current discussion. All that matters is that it satisfies detailed balance, has the Boltzmann 
distribution ip 2 ^ as its stationary state. Also, P is a single-spin flip matrix, i.e. P(S|S') 
vanishes if S and S' differ by more than a single spin. The desired relaxation time of mode 
i is given by 

Ti (L) = -L- 2 /lnA,,(L), (15) 

where Aj is an eigenvalue of Markov matrix P. We obtained the previous expression by 
assuming a single-spin flip Markov matrix, so that the L 2 in the denominator produces a 
relaxation time measured in units of sweeps, i. e. flips per spin. 

IV. TRIAL STATE OPTIMIZATION 

To verify Eq. (|T3|), it is important to obtain estimates that are exact within the range of 
the estimated error. For this purpose we use a set of optimized variational basis functions, 
to which we subsequently apply the projection procedure described in Section 2 to remove 
the variational bias. 

As mentioned, the Monte Carlo projection increases the statistical noise, and the solution 
to this problem is to improve the variational basis functions. We shall now discuss how this 
is done and we consider the problem using the language of the Schrodinger equation. 

We first consider the ground state and review how one can optimize a many-, say 50- 
parameter trial function if)^(R). || The local energy £(R) is defined by 

H^ T (R) = £(R)ip T (R). (16) 

The variance of the local energy is given by 

x 2 = ((H-£) 2 ) = J \MR)\ 2 [£(R)-WdR/ J \MR)\ 2 dR. (17) 

A property that we shall exploit later is that x 2 = for any eigenstate, not just the ground 
state. 

The following sums up the Monte Carlo optimization procedure for a single state: 

1. Sample R\, . . . , R s from ip^ 2 a typical sample size has s ~ 3, 000. 

2. Approximate the integrals in Eq. ( |17D by Monte Carlo sums. 

3. Minimize x 2 as follows, while keeping this sample fixed. For each member of the sample 
Ri, . . . , R s : 
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4. Compute ^t(-Ri), • • • , i^r(Rs)- 

5. Compute TCiPt(Ri), . . . , TCiPt(R s ). 

6. Find S from least-squares fit of 

HMRa) = £MR*), ° 



l,...,s. 



(18) 



7. Minimize the sum of squared residues of Eq. |TH[[] 

This procedure can be generalized immediately to a set of basis functions, as required 
to implement Eq. (|^). The only new ingredient is a guiding function ^ that has sufficient 
overlap with all basis states used in the computation. For this purpose one can conveniently 
use the groundstate raised to some appropriate power less than unity. 

This yields the following algorithm to optimize basis states for n dominant eigenvalues: 

1. Sample Ri, . . . ,R S from 



2. Compute the arrays 



fuW(R 1 )\ 
u^iRt) 



fuW(R s )\ 
u^(R s ) 



J 



(19) 



3. Compute the arrays 



Hu^iRx) 



V 



/ 



[HuV(R s )\ 
HuM(R s ) 



\ 



J 



(20) 



4. Find the matrix elements £u from the appropriately weighted least-squares fit to 



3 = 1 



(21) 



5. Vary the parameters to optimize the fit, as explained below. 

In case of a perfect fit, the eigenvalues of the truncated Hamiltonian matrix E = 
(£y)"- = 1 are the required eigenvalues, but in real life one has to optimize the parameters of 
the basis functions, which can be done as follows: 

1. Divide the sample in blocks and compute one Hamiltonian matrix E per block. 



1 Once the parameters are changed from the values they had in step |], one should use an appro- 
priately weighted sum of squared residues. M 
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2. Minimize the variance of the E-spectra over the blocks. 

The variance vanishes if the basis functions are linear combinations of n eigenstates of TC. 
This gives rise to a computational problem, viz., the variance is near- invariant under linear 
transformation of the This approximate "gauge invariance" gives rise to near-singular, 
non-linear optimization problem. This can be avoided by simultaneously minimizing the 
variance of both the spectrum of the "local" Hamiltonian matrix E and the local energy 8 
of the individual basis functions. 

Finally, the variational bias of the eigenvalue estimates obtained with the optimized basis 
states is reduced by using Monte Carlo to make the substitution discussed previously 

\u {i) ) ^e- HT \u (i) ). (22) 

For this purpose, one has to use the short-time approximation of exp(— rir). To apply 
the preceding scheme to the problem of critical dynamics, all one has to do is to make use 
of the fact the analog of the quantum mechanical evolution is the symmetrized Markov P 
of stochastic dynamics, which is defined as 

£(S|S') = ^yP(S|S')V>B(S'), (23) 

in terms of which we have the correspondence 

e~ HT -> P l . (24) 



V. XE TRIMER: A TEST CASE 

As an example that illustrates the accuracy one can obtain by means of the optimization 
schemes discussed above, we present results for a Xe trimer interacting via a Lennard- Jones 
potential. To be precise, we write the Hamiltonian of this system in reduced units as 



n 



2m 



v 2 + E( 



2)r 



ij i 



(25) 



*<3 



where the r. 
to m 



l3 denote the dimensionless interparticle distances. We define Xe to correspond 
1 = 7.8508 x 1CT 5 , which probably to four significant figures |I(J agrees with Leitner 



et ai. 11 



Table | shows results for variational energies of the lowest five completely symmetric 
states of a Lennard- Jones Xe trimer. The results are compared with results obtained by the 
discrete variable representation truncation-diagonalization method. [11] The basis functions 



used in this computation are of the same general form used in earlier work with an additional 



polynomial prefactor for excited states. [12.13 



Clearly, we obtain consistently lower reduced energies, which we attribute to lack of 
convergence of the results of Leitner et al. |14 
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TABLES 



k 


Ek 


a 


Leitner et al. 





-2.845 241 50 


1 xl0~ 8 


-2.844 


1 


-2.724 955 8 


1 xl0~ 7 


-2.723 


2 


-2.675 065 


1 xl0~ 6 


-2.664 


3 


-2.608 612 


2 xl0~ 6 


-2.604 


4 


-2.592 223 


3 xl0~ 6 


-2.580 



TABLE I. Variational reduced energies compared with estimates of Leitner et al. 



VI. CRITICAL POINT DYNAMICS: RESULTS 

Next we briefly address the issue of the choice of trial functions for the eigenstates of 
symmetrized Markov matrix P. We write 

«(S) = f(S) x ^ B (S). (26) 

For the modes we considered, /(S) was chosen to be a rotationally and translationally 
invariant polynomial in long-wavelength Fourier components of S, the lowest-order one of 
which is simply the magnetization. Corresponding to the order parameter and energy-like 
modes, we considered polynomials either odd or even under the transformation S — > — S. 

We briefly discuss some of the results that illustrate the validity of Eq. QIBp. Figure [I], 
shows plots of the effective amplitudes for the three dominant odd, and two dominant even 
modes of three different Ising models on L x L lattices. Of the three Ising models we 
studied, the first one, the NN model, had nearest-neighbor couplings only. The other two 
also had next-nearest-neighbor couplings. In one of them, the equivalent neighbor or EQN 
model, both couplings were of equal ferromagnetic strengths. In the third or NEQ model, 
the nearest-neighbor coupling was chosen ferromagnetic and of twice the magnitude of the 
antiferromagnetic next-nearest-neighbor coupling. 
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FIG. 1. Universality of relaxation-time amplitudes, shown in a plot of the effective, 
size-dependent amplitudes An on a logarithmic scale. To separate data points for the three mod- 
els, the NEQ data were displaced to the left and the EQN data to the right. The data collapse 
predicted by Eq. ( |i~3|) was produced by fitting the metric factors of the NN and NEQ models. 
Amplitudes of odd and even states alternate in magnitude. 



To obtain estimates of the amplitudes of the relaxation modes, we fit the computed 
correlation times to expressions of the form 

n (L) « L Z Y, a kt L- 2k . (27) 

k = 

In our computation of the non-universal metric factors, this quantity was set equal to 
unity by definition for the EQN model. Table [IT] shows the metric factors computed for each 
mode separately as the ratio of the computed amplitudes. In agreement with the scaling 
prediction in Eq. (|13|), the computed metric factors depend only on the model but not on 
the mode. 

TABLE II. Non- universal metric factors m K , as defined in Eq. (|b|), computed for the NN and 
NEQ models. The modes indicated by ol, o2, and o3 are odd under spin inversion; the remaining 
two, e2 and e3, are even. 



ol 

c2 
o2 
c3 
o3 



NEQ 
2.389(1) 
2.394(2) 
2.393(2) 
2.391(2) 
2.385(4) 



NN 
1.5569 (5) 
1.5569 (5) 
1.5567 (6) 
1.554 (2) 
1.554 (2) 
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Finally we mention that the spectral gaps of the Markov matrix vary over a considerable 
range 

1 - Xi(L) « L^ d+z) « L- 4 - 17 , (28) 

z.e. from approximately 3 x 10~ 3 for L = 4 to 3 x 10~ 6 for L = 21. For details of the 
numerical analysis based on Eq. ( P7|) we refer the interested reader to Ref. ||. Suffice it to 
mention that the value obtained for the universal dynamic critical exponent z featured in 
Eq. ( |T3"D is z = 2.167 ± 0.002 which is indistinguishable from 13/6. 
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